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ABSTRACT 

Results of a 3D MHD simulation of a sunspot with a photospheric size of about 20 Mm are presented. 
The simulation has been carried out with the MURaM code, which includes a realistic equation of 
state with partial ionization and radiative transfer along many ray directions. The largely relaxed 
state of the sunspot shows a division in a central dark umbral region with bright dots and a penumbra 
showing bright filaments of about 2 to 3 Mm length with central dark lanes. By a process similar to 
the formation of umbral dots, the penumbral filaments result from magneto-convection in the form of 
upflow plumes, which become elongated by the presence of an inclined magnetic field: the upflow is 
deflected in the outward direction while the magnetic field is weakened and becomes almost horizontal 
in the upper part of the plume near the level of optical depth unity. A dark lane forms owing to the 
piling up of matter near the cusp-shaped top of the rising plume that leads to an upward bulging of 
the surfaces of constant optical depth. The simulated penumbral structure corresponds well to the 
observationally inferred interlocking-comb structure of the magnetic field with Evershed outflows along 
dark-laned filaments with nearly horizontal magnetic field and overturning perpendicular ('twisting') 
motion, which are embedded in a background of stronger and less inclined field. Photospheric spectral 
lines are formed at the very top and somewhat above the upflow plumes, so that they do not fully 
sense the strong flow as well as the large fleld inclination and signiflcant fleld strength reduction in 
the upper part of the plume structures. 

Subject headings: MHD - convection - radiative transfer - sunspots 



1. INTRODUCTION 

The physical understanding of sunspot structure has 
been hampered for decades by 1) insufficient resolu- 
tion of the fine structure by observations, 2) lack of 
information about the layers below the visible surface, 
and 3) insufficient computational power to perform ab- 
initio 3D MHD simulation of a full sunspot includ- 
ing the surrounding granulation. Recently, we have 
seen considerable progress on all three of these fronts: 

1) adaptive optics, image selection and reconstruction 
at ground-based telescopes and the advent of spectro- 
polarimetry in the visible from space with the Hinode 
satellite have led to a wealth of new information about 
the fine structure of sunspot um brae and penum brae (e.g. 
Bharti et aIll2QQ7l: iLanghans e t al. 2007; Ichim oto et all 
20071 : iRiethmiiller et al.l l2008: Rimmele & Marino "2006") 

2) local helioseismology has star ted to probe the sub - 
surface structure of sunspots (e.g. JCameron et 311120081 ). 
and 3) the ever-increasing computational power of par- 
allel computers have made ab-initio simulations of full 
sunspots come in to reac h. 

While Camer on et al.l (|2007f ) simulated solar pores of 
up to about 3 Mm diameter and did not find indications 
for the development of a penumbral structure, the first 
attempt to simulate a sunspot together with the sur- 
rounding granulation is due to Heine mann et al.l (|2007i ). 
They considered a rectangular section of a (slab-like) 
small sunspot of about 4 Mm diameter. The main result 
of this simulation is the formation of filamentary struc- 
tures in the outer part of the spot, various properties of 
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which (such as dark cores, outfiows, and strongly inclined 
magnetic field) are consistent wi th observational results . 
However, the filaments found by Heinemann et alj (|2007l ) 
are much shorter than the typical lengths of real penum- 
bral filaments and the overall extension of the simulated 
penumbra is very small. 

Here we report about results of a sunspo t simulation 
with the MURaM code (Vogler et al. 2005). The sim- 
ulated sunspot has a total diameter of about 20 Mm, 
shared about equally by umbra and penumbra. 

2. SIMULATION SETUP 

The primary numerical challenges of performing a 
large-scale sunspot simulation including umbra, penum- 
bra and granulation are the significant variations of 
13 — Sttp/^^ and Alfven velocity encountered at a given 
geometrical height. While a very low value of {3 primarily 
imposes stability problems, very high Alfven velocities of 
more that 1000 km-s~^ above the umbra of a spot lead 
to unacceptably small time steps for an explicit code. In 
order to cope with these problems, we have modified the 
MURaM code as outlined in Appendix [Al 

The basic concept of ou r simulation is similar to that 
of lHeinemann et al] (|2007l ): we consider a slender rectan- 
gular section of a sunspot. However, our computational 
domain is considerably wider and deeper than theirs, 
permitting us to simulate a much larger sunspot. The 
rectangular box spans 36.864 Mm x 4.608 Mm in the hor- 
izontal (x, y) directions and 6.144 Mm in the vertical {z) 
direction. The side boundaries are periodic. The top 
boundary is closed for the fiow and the magnetic field is 
matched to a potential field above. At the bottom, the 
boundary is open for th e fiow and the ma gnetic field kept 
vertical (for details, see lVogler et al.|[2QQ 5). In regions of 
strong magnetic field {B > 1000 G), the bottom bound- 
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Fig. 1. — Continuum intensity image at 630 nm of the simulated sunspot and its environment (doubled in the y-direction). The bright 
umbral dots and penumbral filaments have peak intensities between 40% and 90% of the average value outside the spot. The penumbral 
filaments reach lengths of 2-3 Mm. The white frame indicates the filament studied in detail in Sec. 3.2. 
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Fig. 2. — Vertical structure of the simulated sunspot. Shown is a color coding of the square root of the field strength on a cut through 
the simulation in x-z plane dX y = 3.84 Mm. White color corresponds to the maximum field strength of 9kG, black to zero field. The 
whitish line indicates the height of the level reso = 1- The average Wilson depression of the spot umbra is about 450 km. Separated from 
the main sunspot, two pore-like structures have formed at x ~ 31 Mm and x 35 Mm, respectively. 
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Fig. 3. — Large-scale structure of the simulated sunspot: quan- 
tities averaged over the horizontal y-direction as functions of x. 
Upper panel: vertical magnetic field, Bz, at reso =0.1 (solid), 
modulus of the horizontal field component, \Bx\^ at reso =0.1 
(dashed), and normalized continuum intensity at 630 nm (dotted). 
Lower panel: magnetic field inclination with respect to the verti- 
cal at T630 = 0.1 (solid) and horizontal velocity, Vx, at tqsq = 1. 
(dashed). 

ary is closed to avoid outflow instabilities in long runs. 
The radiative transfer is treated in the grey approxima- 
tion. 
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Fig. 4. — Continuum intensity image showing details of penum- 
bral filaments. The vertical lines indicate the positions of the ver- 
tical cuts presented in Fig. [5] 

The simulation was started from a thermally relaxed 
non-magnetic run by introducing a 2D (slab) fleld con- 
flguration in the x-z plane with a width of 5 Mm and a 
strength of 10 kG at the bottom of the box, expanding to 
15 Mm at the top. Using a moderate grid resolution of 
48 km X 48 km x 32 km we run the simulation for about 12 
hours solar time. After a very dynamic adjustment phase 
of about 2 hours, elongated fllaments with dark cores of 
about 2 to 3 Mm length started forming in the periphery 
of the umbra, their heads moving inward. After about 
7.5 hours of evolution we restarted from a snapshot and 
increased the resolution to 32 km x 32 km x 21.33 km. 

3. RESULTS 

For the detailed analysis of the simulation results, we 
consider a snapshot taken 30 solar minutes after the start 
of the high-resolution run. We flrst present results con- 
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Fig. 5. — Vertical cuts in the y-z plane of magnetic field strength 
(left), inclination angle (middle) and horizontal velocity (right) for 
the positions indicated in Fig. |4] (the top row corresponding to the 
leftmost cut). The magnetic field strength is saturated at 4kG, 
the velocity at 2km-s~-'^. The filaments correspond to about 2 Mm 
deep channels of weaker and almost horizontal field and an outfiow 
of a few km-s"-*^. 

cerning the global structure of the simulated spot and 
its average properties and then consider the penumbral 
filaments. 

3.1. Large-scale structure 

Fig. [1] presents a continuum intensity image (A = 
630 nm) of the sunspot and its environment. For better 
visibility, the snapshot has been doubled in ^/-direction. 
The spot umbra has a width of about 10 Mm and is sur- 
rounded on both sides by a penumbra of about 4 — 5 
Mm width, harboring filaments with dark cores. The 
umbra shows umbral dots similar to those found by 
iSchiissler k Voded (|2006l) in a local simulation. 

Fig. [2] displays the square root of the magnetic field 
strength on a vertical cut through the simulation box at 
y = 3.84 Mm. The height expansion of the sunspot due 
to the decreasing gas pressure is well visible. Compari- 
son with Fig. [U shows that the regions corresponding to 
the penumbra (roughly in the ranges x = 9... 14 Mm and 
X = 23... 28 Mm) mostly have a deep underlying magnetic 
structure. In fact, about the same amount of vertical flux 
emerges in the penumbral part of the simulated spot as 
in t he umbra, so that the penumbra cannot be shallow 
(cf . iSolanki fc Schmidt! [1993, ) . 

Fig. [2] also indicates that the magnetic field is rather 
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Fig. 6. — Horizontal and vertical profiles of various physical 
quantities. Top panel: Horizontal cut along the middle line in 
Fig. 0]at height z = 5 Mm, near the r = 1 level in the filament. 
Shown are the vertical field component {Bz, black), the horizontal 
field component {Bx, blue), and the horizontal velocity component 
(vx, green). The vertical field drops to very low values in the 
filaments, resulting in an almost horizontal field of 1 — 1.5kG. Along 
the centers of the filaments, we have an outward directed horizontal 
fiow. Middle panel: Vertical profiles of the same quantities through 
the center of the second filament from the left in the top panel. In 
addition, the red line shows the temperature profile. The horizontal 
outfiow peaks close to the r = 1 level, where the magnetic field has 
the largest inclination. Bottom panel: Vertical profiles at the center 
of the fourth filament from the left in the top panel. 

inhomogeneous between the level reao = 1 (indicated by 
the white line) and about 2 Mm below. This is caused by 
magneto-convection in the form of upflow plumes which 
lead to reduced field strength owing to expansion and flux 
transport by the overturning flow. As we shall discuss 
further below, the underlying processes for the forma- 
tion of penun ibral filaments are the sam e as those for the 
umbral dots ([Schiissler fc VQg"Ierll2006l ). Above reso = 1, 
the field again becomes rather smooth: owing to the low 
plasma P in these layers overlying umbra and penumbra, 
gas pressure gradients cannot maintain strong field inho- 
mogeneities and the field structure has to become almost 
force- free. 

The sunspot has already lost some amount of flux to its 
environment, which shows a plage-like mean vertical flux 
density of about 200 G at tqso = 1- Most of this flux has 
been assembled in small-scale flux concentrations in the 
intergranular lanes (bright structures in Fig.[T]), but also 
pore-like structures have formed, probably supported by 
the periodic boundary condition. 

A more quantitative account of the large-scale struc- 
ture of the simulated sunspot is provided by Fig.[3l which 
shows horizontal profiles of average quantities deter- 
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Fig. 7. — Grey-scale maps of physical quantities on the surface reso = 0.1 for the filament indicated by the rectangular box in Fig. [T] 
The ranges between minimum (black) and maximum (white) of the various quantities are: \Bx\: 670 G ... 1970 G; By-. — 940G ... 640 G; 
Bz-. 980 G ... 3050 G; B inclination with respect to vertical: 17.5 deg ... 61deg; -Vx: -1.4km-s~i ... 3.3km-s~^ Vy: -0.95km-s~i ... 
0.66 km-s"-*^; Vz'. — 2.1km-s~^ ... 1.5km-s~-'^; Iqso- 0.13 ... 1.02 of the average value outside the spot. 



mined in the layers accessible to observation. The profiles 
of the magnetic field components and its inclination with 
respect to the vertic al agree well with the observationa l 
curves determined bv lKeppens fc Martinez Pilletl ([1996), 
although the field strength in the simulated penumbra is 
somewhat high. Average outward horizontal flows with 
peak velocities of about 2km-s~^ are clearly visible in 
both penumbral regions. The main discrepancy between 
the observed and simulated penumbra is the lower mean 
brightness in the simulation (about 60% of the quiet 
Sun) compared to the observed value of about 75% (e.g., 
ISchlichenmaier fc Solankill2QQ3l ). Guided by our experi- 
ence with umbra simulations, we conjecture that this dis- 
crepancy could be due to insufficient thermal relaxation 
of the magnetic region, so that the magneto-convective 
processes providing the structuring and the energy trans- 
port in the penumbra are still not completely developed. 



3.2. Penumbral filaments 

Fig. m shows an enlargement of the penumbral region 
on the right-hand side of the continuum intensity map 
(Fig.[T]). The most conspicuous structures are dark-cored 
filaments of up to a few Mm length. On their end facing 
the umbra, the filaments typically show bright 'heads', 
which propagate into the umbra. The typical lifetime of 
the filaments is about one hour. The three lines in Fig. [4] 
indicate the positions of vertical cuts shown in Fig. [5l 
which give the magnetic field strength, 5, the inclination 
angle with respect to the vertical, arcsin(5a^/|5|), and 
the horizontal flow velocity, Vx. 

It is obvious from Fig. [5] that the penumbral filaments 
correspond to slender structures of about 2 Mm depth 
(near the umbra) and a few hundred km width with 
significantly reduced overall field strength. The verti- 
cal field component, 5^, is reduced to close to zero near 
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Fig. 8. — Profiles of the magnetic field and velocity compo- 
nents at optical depth reso = 0.1 along a cut perpendicular to the 
filament shown in Fig. [7] (cut in y-direction at x = 2 arcsec). [/p- 
•per •panel: Vx (solid), Vy (dashed), Vz (dotted). Lower panel: Bx 
(solid). By (dashed), Bz (dotted). 

the top of the filament, so that an inclination angle of 
almost 90deg results, i.e., the field becomes nearly hor- 
izontal. In the inner penumbra, the filaments develop 
within regions of strong magnetic field; only in the pe- 
riphery of the spot, where the vertical thickness of the 
penumbra drops to less than 2 Mm, the filaments be- 
come more connected to the almost field-free convection 
zone beneath. Therefore, the filaments in this simula- 
tion do not originate from convection penetrating into 
the penumbra from beneath; they are rather similar to 
the plume structures leading; to umb ral dots in the sim- 
ulation of ISchiissler fc Voglerl (|2006l ). here modified by 
the presence of a significant horizontal field component. 

The third column in Fig. [5] shows that the penum- 
bral filaments are associated with strong horizontal flows 
away from the umbra. These flows reach their largest ve- 
locities of a few km-s~^ in the upper and outer parts of 
the filaments, where the magnetic field is most strongly 
inclined. 

Fig. [6] provides a detailed view of the physical condi- 
tions in the filaments by showing horizontal and vertical 
profiles of magnetic field and velocity components for 
the cut shown in the middle row of Fig. [H The horizon- 
tal profiles at z = 5 Mm in the top panel illustrate the 
connection between horizontal outflow and almost hori- 
zontal magnetic field. The vertical cuts at the positions 
y = 1.94 Mm and y = 2.6 Mm, respectively, shown in 
the middle and bottom panels of Fig. [6] indicate that the 



outflows are concentrated in the near-surface layers. It 
is tempti ng; to associate these flows with the Evershed 
effect fcf. IScharmer et al.l[2QQ8l ). 

In what follows, we will study in more detail one 
prototypical filament (outlined by the white square in 
Fig. [T]). Fig. [71 shows a magnified continuum inten- 
sity image at 630 nm (bottom left panel) together with 
maps of various quantities calculated at the surface 
^^630 = 0.1, roughly corresponding to the layer domi- 
nating the spectro-polarimetric information obtainable 
with the often used neutral iron lines at 630.15 nm and 
630.25 nm. This gives a first idea of the actually observ- 
able consequences of the physical processes underlying 
the penumbral filamentation, although definite results 
will require a detailed comparison with synthetic Stokes 
profiles from a non-grey simulation, also taking into ac- 
count image degradation by a realistic point-spread func- 
tion and noise. Such a study is beyond the scope of this 
paper. 

The intensity image in Fig. [71 shows a bright filament 
of a few hundred km width, pervaded by a dark lane of 
about 100 km width. The head of the filament has moved 
some way into the umbra and has nearly disconnected 
from the filament, forming a peripheral umbral dot. The 
maps at constant optical depth tqsq =0.1 show an up- 
fiow {vz > 0) of up to 1.5 km-s~^ at the center of the fila- 
ment, which is connected to downflows at both sides via 
outflows in ±?/-direction, perpendicular to the filament 
axis. In combination with the longitudinal flow along 
the filament, such a flow pattern is consistent with the 
recent observations of ^ twistin g motions' in p enumbral fil- 
aments (jIcliinioto_el_^ The 
strong horizontal flow away from the umbra {—Vx > 0) 
along the dark lane is associated with a weaker, laterally 
more extended, inward return flow at the periphery of 
the filament. The profiles of the three velocity compo- 
nents along a cut in ^/-direction through the filament at 
X = 2 arcsec shown in the upper panel of Fig. [HI clearly 
reveal the flow pattern of a rising plume with laterally 
overturning motion and a strong longitudinal outflow. 

The maps of the magnetic field components in Fig. [71 
and the corresponding profiles along the perpendicular 
cut at X = 2 arcsec (lower panel of Fig. [8|) show a reduc- 
tion of the vertical component, while the component 
along the filament, Bx, stays almost constant. Accord- 
ingly, the field is more strongly inclined and its strength 
reduced in the filament. The component perpendicular 
to the filament axis reverses sign at the filament center 
(above the dark lane), corresponding to the cusp-shaped 
geometry of the top part of the filament: the strong field 
sideways of the filament expands and 'closes' above the 
filament. This configuratio n has recently been confirmed 
observationally by Borrero et al.l ([2008). 

The structure below the visible filament and its rela- 
tionship to the surfaces of constant optical depth that 
are relevant for observations are illustrated by vertical 
cuts perpendicular to and along the filament, which are 
shown in Figs. [9l and [TOl respectively. The perpendicular 
cuts at X = 2 arcsec (Fig. [9]) reveal a structure which is 
very similar t o that underlying; the umb ral dots in the 
simulations of ISchiissler fc Voglerl (|2006f ): a strong up- 
flow plume in the center with narrow downflows at its 
sides leads to an elevation of the iso-r surfaces by about 
200 km, also clearly visible in the temperature profile. 
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Fig. 9. — Vertical cuts at x = 2arcsec perpendicular through the filament shown in Fig. [7] The vertical scale is the same as used in 
Figs. [5] and [G] where zero height corresponds to the bottom of the computational domain. The grey-scale ranges between minimum (black) 
and maximum (white) of the quantities shown are: \B\: 480 G ... 4220 G; B inclination with respect to vertical: 20.1 deg ... 90deg; T: 
3890K ... 15410K; -Vx: LSSkm-s^^ ... 2.09km-s-i; Vy: -O.Tkm-s^^ ... O.Skm-s^^; Vz: -2.3km-s-i ... 2.2km-s-i . The white lines 
indicate the levels of reso = 0.1, and 0.01, respectively. 
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Fig. 10. — Vertical cuts along the center of the filament at constant y = 0.9arcsec in Fig. [71 The grey-scale ranges between minimum 
(black) and maximum (white) of the quantities shown are: \B\: 480 G ... 4220 G; B inclination with respect to vertical: 20.1 deg ... 90deg; 
-Vx-. -1.3km-s-i ... 3.6km-s-^ 1;^: -2.7km-s-i ... 2.5km-s-^ 
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Fig. 11. — Schematic illustration of the basic structure of penum- 
bral filaments as suggested by Zakharov et al. ( 2008) on the basis of 
spectro-polarimetric observations. Shown is a vertical cut perpen- 
dicular to the filament axes. The bright semicircular areas indicate 
the uppermost part of the upfiow plume, the curved arrows repre- 
sent the overturning convective fiow. Circled crosses indicate the 
almost horizontal magnetic field and velocity along the filaments. 
The lines in the grey area are projections of the less inclined mag- 
netic field lines outside the filaments. The intrusion of the plume 
and fiux transport by the overturning motion have pushed aside 
the less inclined field between the filaments. 

The magnetic field is strongly reduced in the plume since 
1) plasma moving upward expands owing to the strong 
pressure stratification and 2) the overturning horizontal 
fiow transports fiux to the downfiow regions at the edges 
of the filament. The latter effect preferentially weakens 
the vertical field component while horizontal field is re- 
plenished by the upfiow. The presence of the strongly in- 
clined magnetic field defiects the central upfiow outward, 
leading to a horizontal fiow of about 2 km-s~^ away from 
the umbra of the spot. 

There is some indication that occasionally part of the 
overturning fiow is recirculated into the upfiow plume, so 
that the fiow pattern becomes reminiscent to roll convec- 
tion as originally suggested by Danielson (1961). To see 
this, consider the velocity components shown in the lower 
panels of Fig.[9l the central upfiow (vz) drives perpendic- 
ular outfiows {vy) near the top of the filament, which feed 
the downfiows adjacent to it. About 300 km deeper, the 
sign of Vy is reversed, so that we now have a perpendic- 
ular infiow, which converges with the central upfiow and 
closes the roll. It is not clear at this moment 1) whether 
this is a robust or a transient feature of the fiow pat- 
tern, 2) whether it is really roll- type convection or more 
related to Kelvin-Helmholtz instability of the downfiow, 
and 3) whether it is realistic at all since the simulation 
still has a much larger effective magnetic diffusivity than 
the plasma in a real sunspot. Anyway, its rather small 
depth extension indicates that such roll convection prob- 
ably does not contribute much to the convective energy 
transport, which is mainly provided by the deep upfiow 
plume. 

Fig. [To] shows physical quantities on a vertical cut 
roughly along the dark lane. It reveals that the structure 
underlying the visible filament is extended in depth as far 
as the upfiows and the reduction of the field strength are 
concerned. On the other hand, the horizontal outfiow 
(vx) and the field inclination sharply peak near optical 
depth unity and lead to the formation of a narrow, al- 
most horizontal fiow channel. 

Altogether, the structure of the simulated penum- 
bral filaments may well be represented by the sketch 
shown in Fig. [TTl which illustrates the essence of 



the high-resolution sp ectro-polarimetric observations of 
IZakharov et"all (|2008[ ). 

A direct connection of the Evershed fiow with the ba- 
sic magneto- co nvective proces s in th e penumbra has been 
suggested by Scharmer et al.l (l20Q8f) o n the basis of the 
simulations of Heine mann et al.l (120071). In fact, we find 
a dominant outward fiow due to the defiection of the up- 
ward fiow in the central part of the filament, which is 
turned outward by the presence of the inclined magnetic 
field. The corresponding return fiow is mainly located in 
the regions with less inclined and stronger magnetic field 
between the filaments (see Figs. [7] and [8]); it therefore 
has a smaller horizontal component at the same (optical) 
depth. Consequently, the observable average horizontal 
velocity is outward (see the lower panel of Fig. [3]) in ac- 
cordance with the observed Evershed effect, even though 
there is no significant overall net outward mass fiux in 
the penumbra. This possibly resolves the long-standing 
problem of the source and disposal of the mass trans - 
ported by the Evershed fiow (e.g., Solanki et al.l 11993 ). 
The existence of the weaker inward return fiows is a pre- 
diction of the simulation that should be testable with 
high-resolution observations. Note that Jchimoto et al.l 
(|2008l ) have found first indications for fiows between the 
Evershed fiow channels. In our current simulation, the 
penumbral filaments appear to be more separated from 
each other than shown by observations. While the re- 
turn fiow is clearly visible around separated individual 
filaments, there are indications that it might be less eas- 
ily detectable in the case of more densely packed fila- 
ments. This needs to be investigated in the future with 
simulations carried out at higher spatial resolution. 

The hot upfiow extends over most of the length of the 
filaments, which provides a rather efiicient energy trans- 
port to supply the radiative losses of the penumbra. In 
contrast, the moving-fiux-tube model with only localized 
upfiows is far too inefiic ient to explain the average bright - 
ness of the penumbra (jSchlichenmaier fc Solankil 120031 ). 
Similarly, the localized plumes underlying umbral dots 
can only maintain the much lower brightness of the um- 
bra in comparison with the penumbra, although the basic 
magneto-convective process is very similar. 

The location of the iso-r lines in Figs. [9] and [10] demon- 
strates that the main part of the upfiow plume and the 
overturning motion underlying the filament is below the 
visible layers. Spectro-polarimetric observations corre- 
sponding to r63o — 0.1 . . . 0.01 just scratch the 'tip of the 
iceberg' and reveal only the uppermost part of the Ever- 
shed fiow and magnetic cusp structure, which is already 
laterally much more homogeneous than the underlying 
main part of the filament. This situation is very simi- 
lar to that of umbral dot observations and in the past 
probably has led, together with the effect of insufficient 
spatial resolution, to ambiguities in the interpretation of 
spectroscopic observations of velocity and magnetic field. 

4. DISCUSSION 

The properties of our simulated sunspot are consistent 
with the general picture of sunspot structure that has 
emerged from observational studies. This applies to the 
overall structure (e.g., distinction between umbra and 
penumbra, average 'radial' profiles of the magnetic field 
components and inclination angle, average outfiow in 
the penumbra) as well as to the detailed properties of 
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the fine structure of umbra and penumbra. The penum- 
bral filamentation results from magneto-convective 
energy transport in the form of hot rising plumes, 
very similar to the proce ss giving rise to umbral dots 
(ISchiissler fc Voglerl l2QQ6'). The inclined magnetic field 
near the periphery of the spot causes a symmetry 
breaking which leads to elongated filaments with strong 
outfiows along fiow tubes of nearly horizontal field 
near optical depth unity. In addition to the fiow along 
the filament, the upfiow also turns over into a motion 
perpendicular to the filament axis. Dark lanes appear 
above the strongest upfiows owing to the upward bulging 
of the surface of optical depth unity and the piling up of 
plasma in a cusp-shaped region at the top of the filament, 
above which the less inclined field outside the filament 
becomes laterally fairly homogeneous. The horizontal 
outfiows are concentrated along the dark lanes. All these 
properties are consistent with recent observational re- 
sults (e .^., [Bellot Rubio e t al. 2005; Rimmele & Marino 
2QQ6I: iL anghans et al.' 2007; Ichimoto et all 120071: 
Borrero et al. 2008,; van Noort fc Rouppe van der VoortI 
20081 : IZakharov et al.ll200a ). 

Our results are also consistent with many properties of 
the sh ort penumbral filaments found bv lHeinemann et al.l 
(|2007l ). who gave interpretations along very similar lines. 
The fact that simulations with two rather different nu- 
merical codes lead to basically the same picture for 
the formation of penumbral structure indicates that, in 
spite of all differences in detail, the simulations have in- 
deed captured essential physical processes. The expla- 
nation for the Evershed effect as a natural consequence 
of ris ing plumes in an inclined field (cf. IScharmer et al.l 
l2008f ) connects this fiow directly to the basic magneto- 
convective structure of the penumbra, so that it should 
occur whenever penumbral structure is present. The ge- 
ometry of the fiow pattern is such that the observable 
average outfiow velocity need not to be connected with 
a net outfiux of mass. This does not exclude the addi- 
tional presence of siphon fiows (e.^., [P egenh a^tl 1 19931 : 
iMontesinos fc Thomasl[l997l : |Solanki et al. 199^ but it 
is much less clear whether the pressure gradients required 
for a sustained outward siphon fiow are maintained al- 
ways and everywhere in all penumbrae. 

What can be said on the basis of our simulation results 
about the various models that have been proposed to ex- 
plain the penumbral structure? First of all, we do not 
see evidence for the 'mo ving flux tube' model and inter- 
change convection (e.g. ISchlichenmaier et al.l 11998a, b). 
In our simulations, the progression of the filament heads 
toward the umbra during their formation phase is not 
caused by the inward motion of a narrow fiux tube, 
but rather due to the expansion of the sheet-like upfiow 
plumes along the filament. 

Furthermore, the simulations co mbine various aspects 
of the 'embedded fiux tube' model (|Solanki fc MontavonI 
119931 : iBorrero et al. 2005, 2006| ) and the 'gappy penum- 
bra' model (Spruit fc Scharmer 2006: Scharmer fc SpruitI 
[2OO6). However, the simulated penumbral filaments are 
neither intrusions of field-free plasma from below nor 
are they confined to almost horizontal fiux tubes dis- 
connected from their environment, particularly in depth. 
The uppermost part of the plume structure forming the 
penumbral filament with its strong horizontal fiow and 
almost horizontal field could possibly be represented by a 



kind of embedded fiux tube. The main feature missing in 
the embedded fiux tube models is the overturning convec- 
tion within the filament and the deep-reaching upfiow of 
plasma that provides the primary energy supply. In this 
respect, the underlying plume structure with its reduced 
(albeit non- vanishing) field strength has much similarity 
with the 'gappy' configuration of the penumbra. This 
scenario captures the convective origin of the penum- 
bral filamentation, even though the gaps form within the 
strong magnetic field. As a consequence, the gaps con- 
tain a horizontal field and, in most cases, are not con- 
nected to the almost field free convecting plasma below 
the penumbra. 

There is no clear evidence in our simulations that the 
penumbral structure is affected or eve n caused by th e 
fiuting instability as suggested by Wei ss et all (|2004l ). 
However, the periodic boundary condition in the hori- 
zontal {x) direction used in the simulation implicitly cor- 
responds to the existence of identical sunspots just about 
20 Mm from the penumbral boundaries, which certainly 
affects the field structure, particularly the inclination, 
in the outer penumbra. Therefore, the simulation possi- 
bly does not well represen t the c onvective pumping effect 
suggested bv lWeiss et al.l (|2004f ) as a mechanism for the 
downward dragging of magnetic fiux in the outer penum- 
bra. This might provide a possible explanation for the 
still rather small extension of the simulated penumbra. 
Observations in fact indicate that penumbrae are often 
suppressed on the side of a sunspot which faces a nearby 
spot of the same polarity. 

Altogether, our results indicate a new level of realism 
in the theoretical modelling of sunspot structure. The 
properties of the simulated penumbral filaments are con- 
sistent with a variety of observational results and pro- 
vide a basis for a physical understanding of umbral and 
penumbral structure in terms of magneto-convective pro- 
cesses. On the other hand, there are still clear discrep- 
ancies between the numerical results and real sunspots, 
so that there is some way ahead to be covered towards a 
completely satisfactory model. Our penumbral structure 
does not yet appear to be fully evolved and the overall 
extension of the penumbra is still somewhat small. The 
average intensity profile indicates that we have simulated 
the development of an inner penumbra, while the outer 
penumbra mig ht be more st rongly affected by convective 
pumpmg (Weis s et all 12004 )). 

The lower boundary condition remains arbitrary since 
we still have no reliable observational constraints con- 
cerning the subsurface structure of sunspots. Computa- 
tional limits have forced us to use a rather coarse spatial 
resolution of 32 km and a fairly small computational box. 
Furthermore, we could only cover a relatively short over- 
all evolution time. As a consequence, the effective diffu- 
sivities in the simulation are still much larger than the 
real values. Test calculations with different resolution 
show that 1) first indications for filamentary structure 
appear already at a horizontal resolution of 96 km and 
2) the reduction of field strength in the plumes increases 
somewhat when we move to a resolution of 24 km. On 
that basis, the fundamental physical process of sheet-like 
plume convection appears to be a robust feature. The re- 
sults will certainly change in detail (and, hopefully, be- 
come even more similar to the observed penumbrae) as 
resolution increases, but we do not expect totally new 
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processes replacing those that we have described here. 

The rapid increase in available computational power 
and the foreseeable progress in local helioseismology will 
soon alleviate some of the limitations of the present ap- 
proach and thus enable us to carry out even more realistic 
simulations. 
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APPENDIX 
NUMERICAL SCHEME 

To alleviate the time step constraint for an explicit code brought about by high Alfven velocity, we limit the strength 
of the Lorentz force in the regions of low /3, such that the resulting Alfven velocity has a given upper bound, Cmax- 

Fl = "^^^ JxB . (Al) 

V ^max + '^A 

Here ^'a denotes the value that the Alfven velocity would have without limitation. The modification of the Lorentz 
force is applied to regions where ^a > Cmax or, in terms of the P value, where 

^^max \ <^max / 

For the simulations presented in this paper, we have used a value of Cgound = 7 km-s~^ to estimate the speed of sound 
above the sunspot umbra and set Cmax = 31 km-s~^. Since the magnetic field is already in an almost force free state 
for /3 < 0.05, this correction has only a minor infiuence on the overall force balance. We have tested this with a 2D run 
for which we varied the value of Cmax by a factor of 3 and found no significant difference. This approach increases the 
explicit time step limit by about 2 orders of magnitude at zero computational expense, in contrast to the alternative 
of an implicit treatment of the Lorentz force. 

The presence of high and low /3 regions as well as high and low Mach number fiows presents a significant challenge 
for a numerical scheme to properly resolve all regimes without being too diffusive in any of them. To this end, we have 
changed the artificial diffusivity scheme of the MURaM code. After a piecewise linear reconstruction of the solution, 
Ui^ where the reconstruction slope in each cell Aui is limited by an appropriate slope limiter, we use the extrapolated 
values at the interface ui = Ui -\- 0.5 Aui and Ur = iXi+i — 0.5 A-u^+i to compute the diffusive fiux 

1 = q+i (/>(^r - Ui.Ui^i - Ui) {Ur - Ui) , (A3) 

where c denotes the characteristic velocity. Using (j) = 1 reduces the scheme to a standard (at least) second order fiux 
(depending on limiter) such as used in a second-order Lax Friedrichs scheme. For the 4th order MHD scheme of the 

MURaM code ^ we found that a choice of = ^ for (ur — ui) • (i^i+i — Ui) > and = otherwise (no artificial 

steepening) represents the best compromise between maximum stability and minimum diffusion. The diffusivity of the 
scheme is controlled through the slope limit er, for which we use a linear combination of the most diffusive Minmod and 
least diffusive Superbee limiter (see, e.g. LeVeque 1990): £ Minmod + (1 — e) Superbee and allow £ to vary as function 
of /3 (typically 0.5 in high- and less than 0.2 in low-/? regions). Furthermore it turned out that it is sufficient to use 
c = 0.1 Csound + ^ + ^aif for the characteristic velocity, which significantly reduces the diffusivity in low Mach number 
fiows. We apply this scheme to all MHD variables and account for effects of (unfortunately unavoidable) mass diffusion 
in the momentum and energy fluxes. The div B error produced by the diffusion scheme is controlled by iterating 

^ = /i(Ax)2grad(divm . (A4) 
ot 

For values of /i = 0.3. ..0.5, typically less than 10 iterations are required to satisfy max {Ax |div5| ) < 10-^5rms, which 
was found to be sufficient. In the results presented in this paper, we did not use any explicit viscosity or magnetic 
diffusivity. The artiflcial diffusivities are added once every time step, after completion of the Runge-Kutta loop of the 
^th Qj^^Qj. MHD scheme. The treatment of diffusivity outlined above leads to a scheme that is fully shock-capturing, 
at least 4th order accurate in smooth regions (higher order is possible depending on slope limiter used), minimally 
diffusive for low Mach number flows, and stable to P values as low as 10~^. 

A very low value of P also leads to problems in codes that use the conservative formulation of the energy equation, 
since the determination of the internal energy requires to compute the small difference between the nearly equal 
values of the total and the magnetic energies. To avoid this problem, we switch to an isothermal equation of state 
in regions where £^int < 10~^£^mag and prevent too small p values by imposing a lower limit for the density when 
Eint < 10~^£^mag- An alternative to this procedure would have been to directly solve an equation for the internal 
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energy at the expense of loosing the advantage of the conservative formulation. Since most of the dynamics are driven 
beneath the photosphere and stability problems only occur in layers with an average density at least 3 orders of 
magnitude less than the photosphere (and therefore of little dynamical importance), we have preferred to stay with 
the conservative formulation. 
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